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Model-independent parametrisations for examining departures from General Relativity have been 
increasingly studied over the past few years. Various observables have been used to constrain the 
parameters and forecasts for future surveys have been carried out. In one such forecast [I], galaxy 
cluster counts were used to constrain the parameters. Here, we carry out a limited set of A-body 
simulations, with a modified Poisson equation, to examine the accuracy of existing mass functions 
for modified gravity cosmologies. As well as altering the gravitational calculation, we include the 
effect of a screening scale to ensure consistency of the theory with solar system tests. Our results 
suggest that if a screening scale exists its effect can be taken into account in the cluster count 
calculation through its effect on the linear matter power spectrum. If this is done, the accuracy of 
the standard mass function formalism in modified gravity theories with reasonably small departures 
from General Relativity, as tested in this work, is comparable to the standard case. 

PACS numbers: 04.80.Cc,98.65.Cw,98.80.-k 



I. INTRODUCTION 

Several explanations have been put forward to explain 
the apparent accelerating cosmological expansion. The 
most popular idea is that the acceleration is due to a dis- 
tinct form of energy in the universe, generically known 
as dark energy. The simplest example of what could 
constitute dark energy is a cosmological constant A, a 
non-dynamical form of energy that may be interpreted 
as arising as an unconstrained integration constant in the 
equations of motion or from the vacuum energy generated 
by the quantised fields that make up the matter content 
of the universe. Dynamical forms of dark energy involv- 
ing scalar fields that are evolving with time can also lead 
to late-time acceleration. Some have also argued that 
the observation that the expansion is accelerating results 
from the misuse of the Cosmological Principle in inter- 
preting the expansion history from the inhomogeneous 
distribution of matter around us. An alternative perspec- 
tive to the above both lines of thought above is to explain 
the departure from the expected late-time behaviour by 
introducing modifications to the theory of General Rela- 
tivity (GR) itself. 

There has been an increased interest in modified gravity 
theories in cosmology over the past decade or so. The 
large number of proposed alternatives to General Rel- 
ativity (GR) has led to interest in examining possible 
deviations from GR in a more model independent way. 
Much work has gone into creating these parametrisations, 
examining their consistency, constraining them with cur- 
rent data and forecasting future experiments (see [5] for 
a comprehensive review). Many different observations 
and combinations of observations have been considered 
for constraining these parametrisations. The key point 
with the combinations is the issue of distinct degenera- 
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cies in different sets of observations. In particular, obser- 
vations such as the integrated Sachs- Wolfe effect in the 
Cosmic Microwave Background (CMB) and cosmic shear 
surveys, both being effects on the energies and trajecto- 
ries of photons that are relativistic, constrain the sum of 
the two scalar potentials in the metric 'J + <E> and thus 
constrain a particular combination of the Modified Grav- 
ity Parameters (MGPs) as will be introduced later. This 
degeneracy is best broken with an observation that de- 
pends on a different linear combination of the scalar po- 
tentials. Observations relating to the growth of structure 
in non-relativistic matter, such as redshift space distor- 
tions or galaxy surveys, are a natural choice for breaking 
the degeneracy. Recently pQ, we argued that combining 
future cluster counts with CMB and cosmic shear ob- 
servations could yield high precision constraints on the 
simplest form of MGPs particularly in the case where 
the modifications are constant at low redshift. 

Galaxy cluster counts can be predicted for a particular 
combination of cosmological parameters by starting from 
the linear matter power spectrum and using mass func- 
tions to calculate how many bound structure will form 
at different mass scales. The mass functions are semi 
analytic formulas whose form derives from the original 
Press-Schechter ideas [3] and its extensions [3HZ|- The 
functions are now calibrated using detailed and exten- 
sive Af-body simulations. In this paper, we present the 
results of a limited set of of Af-body simulations aimed at 
investigating the accuracy of existing mass functions for 
the MGP formalism. Due to the computational demand 
of running the required A^-body simulations we have re- 
stricted our investigation to a single point in the phase 
space of MGPs i.e. a single choice in the modification of 
the strength of gravity that is constant at low redshift. 
As such we are not attempting to accurately calibrate an 
extension of the mass function formalism in the phase 
space of MGPs. We are instead addressing the question 
of whether the existing mass function formalism can be 
used as is in forecasts involving cluster counts such as 
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that carried out in pQ. 

This paper is organised as follows. In section [IT] we dis- 
cuss the mass function formalism and some of the results 
in the literature regarding modified gravity models. The 
modifications required to run A-body simulation in such 
models of parametrised modified gravity are summarised 
in section |III| and the results of running our suite of sim- 
ulations are described in section ITVl We conclude with a 
discussion of the results in section Ivl 



II. MASS FUNCTIONS 

Galaxy clusters are some of the largest collapsed struc- 
tures in the universe. According to the standard scenario 
of hierarchical structure formation within the ACDM 
cosmological model, clusters typically consist of hot gas 
bound in a large cold dark matter halo. Clusters are 
a useful cosmological probe as their size corresponds to 
scales near the linear to non-linear transition in the un- 
derlying dark matter power spectrum. This has several 
consequences: they probe the tail of the matter pertur- 
bation spectrum and are therefore a sensitive probe of 
growth. In addition, galaxy cluster counts can be pre- 
dicted from linear theory, using semi-analytic formulae 
or formulae calibrated from A-body simulations. Cluster 
counts and the properties of clusters have been used to 
constrain cosmological parameters and in particular they 
have been used to constrain dark energy, which may be 
responsible for late-time acceleration of the cosmological 
expansion (see for example |8,'9j), and some studies look- 
ing at constraining the growth factor 7 (see for example 

IE]). 

The number of clusters observable over a fraction of the 
sky / s ky and with a redshift dependent mass resolution 
limit Afii m (z) in a redshift bin spanning the interval z 
to z + Az can be calculated by integrating the comoving 
number density dn/dM of objects with mass M 
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where dV/dzdQ, = r 2 (z)/H(z) is the comoving volume at 
redshift z in a flat universe, with H(z) the Hubble rate 
and r(z) = J Q dz' /H(z') the comoving distance to that 
redshift. 

The comoving number density of objects in a given mass 
range dn/dM is known as the mass function and much 
work has gone into predicting its shape for a given linear 
power spectrum. Early, semi-analytical estimates of the 
mass function resulted in the Press-Schechter formalism 
[3], where the mass function is given by 
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where p is the background density of dark matter today, 
S c = 1.686 is the critical density contrast for collapse 



of a spherical perturbation and o 2 M is the variance of 
the dark matter fluctuations in spheres of radius R — 
(3M/4-7T/9) 1 / 3 defined by the integral of the linear matter 
power spectrum P(k) over wavenumber k 
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with top-hat filter function 



Successive studies have refined the formalism resulting 
in more complex expressions with increased accuracy [3J- 
7 . The functional forms of these mass functions have 
been used to fit the results of A-body simulations for 
ACDM cosmologies and agreement is at the 5-10% level 
for ACDM cosmologies [T2l IT3"] . The inaccuracy increases 
to 10% or more for more general cosmologies [T^l H5] . In 
the what follows we select the following form from [14] . 
henceforth W05: 
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The mass function formalism has been used extensively in 
constraining extensions to the standard ACDM cosmol- 
ogy. Some examples are redshift dependent dark energy 
equation of state models wCT>M [T21[T3] and some models 
of modified gravity including f{R) [HJCEB], DGP QH HE] 
and coupled scalar field cosmologies [THIHO]. However, 
in some cases it is unclear how well the mass function 
formalism will reproduce the correct cluster counts if the 
modifications alter significantly the process of structure 
formation. 

There are three reasons why the expressions commonly 
used to predict cluster counts may not be valid in mod- 
ified models. Firstly the linear matter power spectrum 
the formalism is based on will, in general, differ from the 
standard cosmology for the same matter content. Sec- 
ondly the expressions are calibrated with respect to stan- 
dard ACDM A-body simulations or semi-analytical ar- 
guments and the coefficients obtained in this way may 
not be the correct one in modified cosmologies. As such, 
even if two models may give the same linear matter power 
spectrum, they could differ in their cluster counts predic- 
tions because, for example, the critical density S c may be 
different in the modified model. Thirdly, the expressions 
themselves may not be a corrected parametrisation of the 
transformation of a linear power spectrum into a descrip- 
tion of the mass distribution of clusters. 
Typically only the first possibility is addressed in works 
where cluster counts are used to constrain modified cos- 
mologies. This is achieved by inputing into the formalism 
a modified linear matter power spectrum that has been 
obtained by solving the linear growth equations for the 
dark matter. This is by far the simplest way to account 
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for the modifications since the second and third concerns 
require the use of ,/V-body runs to be addressed. 
Here we will focus on how well the mass formalism does in 
reproducing parametrised modified gravity cluster counts 
if the expressions are unchanged and the modifications 
are solely taken into account by modifications of the lin- 
ear matter power spectrum. The results will inform us 
on whether using the standard formalism can be an in- 
formative tool in investigating the constraining power of 
cluster counts in modified gravity models. In particular 
we will be able to determine whether forecasts made in 
PQ may be biased because of the way the mass function 
formalism was used in the prediction of cluster counts in 
modified gravity models. 

III. N-BODY SIMULATIONS AND MODIFIED 
GRAVITY PARAMETRISATION 

Numerous sets of phenomenological "Modified Gravity 
Parameters" (MGPs) have been suggested in the liter- 
ature, see e.g. [3T] for a partial translation table and 
[22] for a discussion of the differences with some of the 
parametrisations. Most of the parametrisations are phe- 
nomenological modifications to the Einstein equations 
and typically involve a parameter relating to the strength 
of gravity and a parameter relating the two scalar poten- 
tials in the metric. Our potentials are defined as scalar 
perturbations of a flat, FRW metric 

5 oo = -[l + 2¥(f,f)] , 

gij = a 2 (t)[l-2<S>{x,t)}5 lJ , (6) 

where we have made the conformal Newtonian gauge 
choice to fix the remaining two scalar degrees of free- 
dom in the perturbed metric. In ([(jf, is the Newto- 
nian potential and is responsible for the acceleration of 
massive particles whereas $ is the curvature potential, 
which also contributes to the acceleration of relativistic 
particles. For this work, since iV-body simulations are es- 
sentially Newtonian plus an expanding background, we 
have only used one MGP, p. This is the parameter in the 
Poisson equation that controls the strength of gravity. In 
fourier space, 

fc 2 *(a, fc) = ~4:TTGa 2 fi(a, k)p(a)A(a, k) . (7) 

where, a is the FRW scale factor, G is Newton's constant, 
p is the background density of cold dark matter and A 
is the gauge invariant density contrast given by 



where the cold dark matter density contrast is defined as 
S = Sp(a, k)/p(a), v is velocity of the dark matter and H 
is the Hubble parameter. 

In principle, p could be any function of time and space, 
however here we have kept a simple treatment of it. In 



particular, we are considering the situation where, at a 
set redshift of z m g, the parameter transitions from its 
GR value of 1 to some other value. Except for the in- 
clusion of a screening scale, there is no scale dependence 
inserted into the phenomenology. Screening mechanisms 
are an important component of modified gravity mod- 
els as they allow models that differ from GR to repro- 
duce solar system tests. There are several mechanisms in 
the literature, including the chameleon mechanism [23] 
that operates in f(R) gravity and the Vainshtein mech- 
anism 24 that operates in DGP gravity. To further our 
phenomenological approach to modified gravity, we use 
a simpler treatment where the modifications to gravity 
cease below a set length scale, £ sc . 

In order to investigate the use of mass functions when 
forecasting cluster counts in models of parametrised mod- 
ified gravity, we have made some simple modifications to 
the publicly available iV-body, Smooth Particle Hydrody- 
namics (SPH) code Gadget-2 [25]. Gadget-2 is a parallel 
TV-body code that calculates the gravitational accelera- 
tion with a TreePM approach [5S] . For this investigation 
we are only interested in the dynamics of the dominant 
dark matter and we do not include any component with 
non-zero pressure. Since the ./V-body simulation probes 
the low-velocity, Newtonian regime on sub-horizon scales, 
the only modification of gravity that affects the simula- 
tions is any modification to the Poisson equation that en- 
compasses all gravitational aspects of the interaction in 
such codes. In Gadget-2 the Poisson equation is solved 
differently on large and small scales. On large scales, 
a regular grid (the "mesh") is placed in the simulation 
volume. The mass of each particle is then interpolated 
onto the mesh. This discrete density field is then Fourier 
transformed and the Poisson equation is solved in Fourier 
space. This yields the gravitational acceleration of par- 
ticles in the simulation due to distant mass distribution 
in the simulation volume. On smaller scales, the force on 
each particle is obtained by calculating the Newtonian 
potential directly. A hierarchy of increasingly higher res- 
olution cells (the "tree") is constructed. Far from the 
point where the potential is being calculated, it is suf- 
ficient to use cells of coarse resolution containing mul- 
tiple particles treated together at their centre of mass. 
Nearer to the point where the potential is being calcu- 
lated, increased resolution is required until, at sufficiently 
short distances, the effect of each particle is included sep- 
arately. The two components of the potential are com- 
bined via a matching filter to obtain the potential govern- 
ing the overall force acting on each particle at every time 
step. Thus modifications to gravity must be encoded 
in both contributions to the total potential i.e. in both 
Fourier domain solution and direct force calculation. 

Our modifications to Gadget-2 are as follows: Both the 
potential calculated by the particle mesh and the force 
calculated by the tree structure are modified by a fac- 
tor of p for distances above a screening scale of proper 
length £ sc . This screening scale can be switched off such 
that the only modification acting is the change in the 
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TABLE I: Parameters for simulations. 



strength of gravity. We ran the modified gravity sim- 
ulations both with and without the screening scale for 
a series of simulation boxes each with the same volume 
and mass resolution (see table [I] for the parameter values 
used in the simulations). The expansion history in the 
simulations is that of a standard ACDM cosmology for 
both GR and MGP runs. This is in order to isolate the 
effect of the MGP on the growth of structure from any 
change due to the expansion history. Each simulation is 
started at a redshift z; = 50 and average the result of 
each simulation over a number of runs with independent 
random seeds to reduce the effect of sample variance in 
our final cluster counts. 

The linear matter power spectra today used for the mass 
function prediction for both the GR+ACDM runs and 
the modified gravity runs were computed using the mod- 
ified version of CAMB [37| MGCAMB [25] . We have 
modified MGCAMB to incorporate the same modifica- 
tions to gravity as the TV-body simulation, namely a con- 
stant value for (j,, different to the fiducial GR value of 
1, that switches on after redshift z mg and with the op- 
tion to recover GR below a screening scale £ ac . For all 
simulations in this work we have chosen a fiducial value 
of fi = 1.05 for the modification. The second MGP, 
7] = < E > /5', has been kept at its GR value in the MG- 
CAMB code and does not come in to the Af-body simu- 
lations. We consider a flat cosmology with the following 
parameters. The dimensionless Hubble rate in units of 
100 Km s -1 Mpc -1 , h = 0.71, the density of matter 
(baryons + dark matter) and dark energy in units of the 
critical energy density, TO = 0.265, and tt\ — 0.735 re- 
spectively, the optical depth to recombination r = 0.088, 
the amplitude of primordial, super horizon curvature per- 
turbations A s = 2.157 x 10" 9 at k = 0.05 h Mpc -1 and 
their spectral index n s = 0.963. These parameters cor- 
respond to the WMAP 7-year best-fit parameters [2"9"j . 
This model yields a large scale structure normalisation 
of eg = 0.804 for the standard deviations of fluctuations 
on scales of 8 ft -1 Mpc. 

The initial condition for each of the A-body simula- 
tions were calculated using the transfer functions ob- 
tained from the MGCAMB run. The transfer functions 
are used to generate a grid of particles with velocities 
and displacements consistent with the power spectrum of 
the matter expected at the starting redshift Z\. For this 
work we have used the 2LPT [301 131] package to obtain 
the initial velocities and displacements. The package uses 
a second order Lagrange perturbation scheme to evolve 
the displacements and velocities of particles to the de- 
sired redshift. This is more accurate than using the first 
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FIG. 1: The number of clusters as a function of mass of the 
cluster found in the our simulations for the GR case (/i = 1) 
compared to the modified gravity case, with and without 
screening scale. In the modified gravity case the 5% increase 
in /x leads to an increase in clusters on all scales when a screen- 
ing scale is not included. The result is more complicated when 
a screening scale is included with less cluster being produced 
at the low mass scales as expected. 



order Zeldovich approximation. 

We obtained independent initial conditions for each of the 
several realisations run for each box and particle num- 
ber. In the simulations used here, all of the particles are 
cold dark matter particles. For each box size and particle 
number, we also specified the smoothing scale e. The sim- 
ulations all had a starting redshift of = 50, however, 
we ran several simulations starting at higher redshifts to 
check there was no effect on the final results. We have 
also verified that our ACDM simulations reproduced the 
results of earlier papers, in particular the results from 

[3H33J. 

Once the simulations are run the challenge is to identify 
clusters and their masses in order to determine the mass 
function of clusters for the simulated volume. There are 
two main algorithms to find dark matter halos in cosmo- 
logical A^-body simulations, namely spherical overdensity 
algorithms {Mj and friends-of-friends algorithms [351 13S] • 
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FIG. 2: The cluster count found in the GR compared to mod- 
ified gravity and modified gravity with a screening scale of I 
Mpc for the simulation volume. We have divided out the 
respective prediction for GR and modified gravity based on 
the mass function formalism applied to the respective linear 
power spectra. The input spectrum is the same for both mod- 
ified gravity cases in this figure. The fact that the prediction 
does not take into account the screening scale leads to an 
over-prediction of clusters on all scales. The solid line gives 
an indication of the expected sample scatter for the simula- 
tion based on the GR counts and is included as a rough guide 
to the significance of any discrepancy. 



FIG. 3: The ratio of cluster counts found in the modified grav- 
ity simulations with respect to the mass function prediction. 
Both counts are for the same set of simulations that included 
a screening scale of 1 Mpc. The green (crosses) points do not 
include the effect of the screening cut in the linear matter 
power spectrum used in the mass function prediction whilst 
the red (triangles) do. The result shows that the mass func- 
tion formalism gives an accurate prediction of the counts if 
the screening scale is taken into account in the prediction. 



IV. RESULTS 



We used a friends-of-friends halo finder 1 with the stan- 
dard linking length equal to 20% of the inter-particle sep- 
aration, ibox/ v/iVpart. In principle, this linking length 
should be changed to take into account how virialisa- 
tion is different for different cosmologies and gravity the- 
ories, however it is not clear how this should be done 
for parametrised modified gravity and we have made the 
standard choice in all of our calculations. It is now known 
that a systematic bias arises when determining the mass 
of halos with the friends-of-friends algorithm [2]. As 
shown in [14] this can be corrected for by modifying the 
number of particles in each halo according to the pre- 
scription N COII = N X (1 — iV~ - 6 ), where N is the orig- 
inal number of clusters found by the friends-of-friends 
algorithm in each mass interval and N COXI is the unbiased 
value. We have applied this correction to the halos found 
in our simulations and also imposed a minimum number 
of (uncorrected) particles per halo of 40. The halos were 
then binned into logarithmic mass bins of width 0.2, be- 
tween 1O 13 /i _1 M and 10 15 h~ 1 M & where M represent 
a solar mass. 
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In this section we examine the output of the simulations 
at a z = 0. Firstly, we examine the change to cluster 
counts due to the modification of gravity. Figure [l] shows 
the total number of clusters found in the GR simulation, 
as well as the modified gravity cases with and without a 
screening scale. We plot the average value of the 10 inde- 
pendent realisations. The error bars show the standard 
deviation estimated from the 10 samples. As expected, 
the number of clusters increases due to the gravity be- 
ing stronger in our fiducial modified model with fj, > 1. 
The increase is larger towards the higher masses. The 
screening scale reduces the increase in clusters in mod- 
ified gravity as we would naively expect. In addition, 
there is a slight decrease in the number of cluster at the 
lowest masses. It is not clear if this decrement is signifi- 
cant or physically relevant but it may be consistent with 
the fact that at some mass scale the effect of screening 
will mean that objects that would have formed clusters in 
the GR case end up bound together with larger clusters 
which grow faster on scales beyond the screening scale. 
In figure [2] we compare the numbers of clusters in our 
simulations to the predictions in each case. The linear 
matter power spectra used as inputs to the W05 mass 
function are computed using the MGCAMB code, with 
l-i = l for the GR case, and [i = 1.05 for the modified 
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gravity case. For the modified gravity simulations with 
and without the screening scale, the linear mass function 
was calculated without the screening scale, as would be 
the case for models of modified gravity where the screen- 
ing scale is an entirely non-linear phenomenon. The mass 
function works as well for the unscreened modified grav- 
ity case as it does for GR. However, for the simulations 
that incorporate a screening scale, the mass function pre- 
diction no longer works. By inserting the screening scale 
into calculation of the linear matter power spectrum we 
can reconcile the mass function formalism with number 
counts observed in the simulations. The effect of taking 
the screening scale into account in the linear power spec- 
trum used as input into the mass function is shown in 
figure [3j 

V. DISCUSSION 

We have shown that the mass function formalism can 
work for simple phenomenological models of modified 
gravity without a screening scale. To model modified 
gravity with a screening scale, its effect needs to be in- 
cluded in the calculation of the linear matter power spec- 
trum in order to retain the accuracy of the mass func- 



tion formalism. Of course, the mass function formalism 
is likely to break down for more extreme departures from 
GR, however it isn't obvious that large deviations are 
allowed by current data. In addition, more complicated 
modelling of the \x parameter, such as including scale de- 
pendence may reduce the accuracy of the mass function 
predictions. However, as long as the scale dependence 
manifests in the linear matter power spectrum and the 
deviations from GR are relatively small then the mass 
function prediction may still be accurate. This will re- 
quire further work. The work presented here needs to be 
extended in other ways as well, to more precisely deter- 
mine the range of validity of the mass function formalism. 
Varying the cosmological parameters, linking length and 
parametrisation of /i will further test the validity of the 
mass functions. In addition, running bigger boxes and 
boxes with improved resolution will test the mass func- 
tion over a larger range of cluster masses. The mass func- 
tion should also be tested over a range of redshifts and 
with baryons included, rather than solely dark matter 
particles. In addition, an examination of spherical col- 
lapse in model independent modified gravity and/or an 
examination of virialisation in these theories may shed 
further light on finding halos and how mass functions 
work in modified gravity simulations. 
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